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Abstract  A number  of  approaches  to  the  problem  of  using  adaptive  mesh  refine- 
ment in  large  eddy  simulations  are  considered  and  tested.  These  include 
unstructured  and  structured  adaptive  grids,  various  treatments  of  the 
convective  terms  and  a number  of  flow  diagnostic  procedures.  The  ap- 
proaches are  exemplified  on  a rotating  channel  flow. 


Introduction 

An  important  aspect  in  applying  large  eddy  simulation  (LES)  to  tech- 
nologically interesting  flows  is  the  control  of  the  computational  costs 
involved.  Complex  flows  usually  present  a variety  of  regions  in  which 
different  resolutions  are  called  for.  Moreover,  these  regions  are  not  static 
in  the  course  of  the  flow  evolution.  Adaptive  mesh  refinement  is  a nat- 
ural candidate  for  keeping  computational  costs  down  in  these  types  of 
applications.  In  this  paper  a number  of  approaches  are  compared  in 
an  effort  to  provide  information  on  the  benefits  of  applying  adaptive 
mesh  refinement  techniques  to  LES  computations.  A first  comparison  is 
made  between  an  unstructured  grid  algorithm  and  a structured  grid  al- 
gorithm. The  same  computational  procedures  are  implemented  in  both 
codes  and  comparable  resolutions  are  achieved  on  a moderately  complex 
flow.  The  results  are  compared  to  one  another  and  a direct  numerical 
simulation  (DNS)  to  permit  evaluation  of  grid  quality  effects.  An  ad- 
ditional comparison  is  made  between  various  criteria  of  controlling  the 
mesh  adaptation.  A variety  of  flow  feature  diagnostic  procedures  are 
allowed  to  control  placement  of  additional  grid  points.  These  include 
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comparison  to  predictor  steps  on  coarser  grids  and  refinement  according 
to  local  values  of  vorticity. 

1.  Adaptive  Mesh  Refinement  for  LES 

1.1  General  observations 

Adaptive  mesh  refinement  (AMR)  [5]  has  shown  itself  to  be  a powerful 
technique  in  the  computation  of  flows  with  dynamic  localized  structures 
that  require  much  greater  resolution  than  that  needed  for  most  areas 
of  the  flow.  AMR  has  been  widely  applied  to  flows  involving  shocks 
or  chemical  reactions.  There  has  been  less  work  in  the  area  of  applying 
AMR  to  turbulent  flow.  AMR  is  not  suitable  for  turbulent  flows  in  which 
fine  resolution  is  required  everywhere.  There  are  however  a number  of 
flows  in  which  small-scale  structures  appear  and  then  drive  the  overall 
flow.  Such  structures  have  limited  spatial  and  temporal  extent.  The 
streaks  that  appear  in  the  boundary  layer  of  channel  flow  are  an  example. 
The  overall  accuracy  of  the  flow  computation  is  sacrificed  if  the  grid 
resolution  is  not  fine  enough  to  capture  the  streaks.  Extending  the  fine 
grid  required  to  resolve  wall  streaks  to  the  entire  domain  is  prohibitively 
expensive.  The  situation  bears  some  resemblance  to  those  in  which  AMR 
has  been  applied  with  success. 

There  is  however  a significant  difference  between  standard  AMR  appli- 
cations and  turbulent  flow.  In  standard  AMR  the  computational  grids 
are  set  up  quite  frequently,  every  few  time  steps.  The  standard  tech- 
nique used  to  identify  regions  in  which  additional  resolution  is  required 
is  to  compare  predictions  on  a coarser  and  finer  grid  and  use  Richardson 
extrapolation  to  estimate  whether  the  local  truncation  error  is  within 
acceptable  limits.  This  implicitly  assumes  that  the  finest  grids  in  the 
computation  can  achieve  enough  resolution  so  that  the  asymptotic  er- 
ror estimates  based  upon  the  formal  order  of  accuracy  of  the  numerical 
method  are  valid.  In  the  context  of  turbulence  computations  this  is 
equivalent  to  stating  that  the  finest  grids  can  achieve  DNS-like  reso- 
lution. This  is  prohibitive  and  a modification  of  AMR  to  account  for 
unresolved  scales  of  motion  is  necessary. 

Previous  work  has  also  sought  to  address  the  problem  of  differentiat- 
ing between  scales  that  should  be  resolved  and  others  that  can  be  eco- 
nomically modeled.  An  early  idea  due  to  Voke  [18]  was  to  use  multiple 
meshes  as  a tool  in  accelerating  convergence  to  a quasi-steady  state.  The 
Dynamic  Multilevel  (DML)  method  of  Dubois,  Jauberteau  and  Temam 
[7]  can  be  seen  as  an  alternative  to  LES.  Instead  of  modeling  subgrid 
stresses  using  a physical  model,  in  DML  the  small  scales  are  computed 
less  accurately  using  lower-order  schemes  with  larger  time  steps.  Terra- 
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col,  Sagaut  and  Basdevant  [16]  have  suggested  freezing  the  small  scales 
while  a number  of  time  steps  are  taken  on  the  coarser  grids.  In  order  not 
to  lose  small-scale  dynamics,  they  suggest  a periodic  prolongation  of  the 
large-scale  velocity  fields  using  small-scales  stored  from  previous  time 
steps.  Both  the  large  and  small  scales  are  advanced  in  time  for  a few 
steps  after  which  a new  large  scale  velocity  field  is  obtained  by  filtering. 
This  procedure  raises  the  question  of  the  equilibrium  between  small  and 
large  scales  during  the  prolongation  stage.  This  step  of  the  algorithm 
might  effectively  act  as  an  external,  non-physical  excitation.  Quemere, 
Saguat  and  Couaillier  [15]  presented  a procedure  in  which  patched  grids 
of  different  resolutions  were  used  in  a channel  flow.  They  show  that 
considerable  economy  may  be  obtained  by  using  fine  grids  only  in  the 
wall  region  but  also  remark  that  numerical  artefacts  arose  at  the  grid 
boundaries. 


1.2  Governing  equations 

The  governing  equations  are  the  filtered,  compressible  Navier-Stokes 
equations  for  an  ideal  gas 

90  m m dh  = 

dt  Ox i 0x2  dx 3 

U = [ p pui  pu2  pu3  pe  ]J  , 
put 

pupil  + pSn  - Til  - pSji 
pu,U2  + P$i2  - Ti2  - pSj2 
puppi  + p8i3  - ns  - pSi3 

{pe  + p)v,i  -Qi-  pSijUj  - k 


F,  = 


09 

dxi 


with  " denoting  grid  filtering,  and  ~ density- weighted  averaging.  The 
filtered  equation  of  state  is  p — pRT , the  diagonal  term  of  the  subgrid- 
stress  tensor  is  neglected  [8],  the  resolved  energy  is 

pe  = pcvT  + ip  {u\  + u\  + «!)  . 


the  Sutherland  relation  for  constant  Prandtl  Pr  = cpp{0)/k{6)  = 0.7  is 
applied,  and  the  system  is  closed  with  variable  density  eddy-viscosity 
and  diffusivity  models 


Tij  — pVtSij, 


Qi  = 


_ vt  do 

P Pit  dxi  ' 
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1.3  Filtering  procedures 

The  general  multi-level  filtering  framework  proposed  by  Germano  [9] 
is  adopted.  Each  grid  level  has  an  associated  filter  operation  Gn.  Filters 
from  the  same  family  (i.e.  top  hat  in  the  applications  presented  here)  are 
used  at  all  levels.  Filtered  variables  may  be  defined  at  each  grid  level  by 
(<f>) n (x,t)  = ( Gn  * <f>)  (x,  t).  The  coarsest  grid  level  is  denoted  by  0 and 
the  finest  existent  at  a given  space-time  locale  by  L.  The  Navier-Stokes 
equations  from  above,  dfU  -1-  N(U)  = 0 become 

dt(U)n  + N((U)n)  = -Tn 

after  filtering.  The  closure  procedure  parallels  that  presented  for  a single 
filtering  operation  in  the  previous  section. 

1.4  Flow  feature  diagnostics 

Standard  AMR  is  typically  viewed  as  a black  box  that  may  applied  to 
the  solution  of  any  time  varying  problem  that  exhibits  localized  features. 
The  criterion  governing  grid  refinement  is  mathematical  in  nature  and 
typically  does  not  include  any  physical  knowledge  about  the  particular 
problem  being  solved.  On  the  other  hand,  most  multi-level  techniques 
that  have  been  proposed  for  turbulence  simulation  rely  heavily  upon 
physical  modeling  of  the  subgrid  scale  effects.  For  instance,  in  DML  [7] 
the  freezing  of  the  small  scales  when  taking  coarse  grid  time  steps  is 
justified  by  the  different  characteristic  times  in  which  small  and  larger 
scale  turbulence  achieves  local  equilibrium.  The  point  of  view  espoused 
in  this  paper  is  that  probably  both  ingredients  are  required  in  order 
to  achieve  a successful  algorithm.  The  particular  method  studied  here 
is  a combination  of  a posteriori  error  analysis  combined  with  physical 
analysis. 

Initially  the  standard  technique  of  error  estimation  based  upon  Richard- 
son extrapolation  [5]  was  tried.  This  was  unsatisfactory  as  it  led  to  a 
rapid  increase  in  the  overall  number  of  points.  A resolution  typical  of 
DNS  was  set  up  by  this  procedure.  The  approach  that  was  successful 
consists  of  a combination  of  error  estimation  and  physical  reasoning. 
Recently,  Adjerid  et  al.  [1]  have  analyzed  the  error  of  a class  of  discon- 
tinuous Galerkin  methods  applied  to  hyperbolic  problems  that  include 
the  standard  finite  volume  schemes  typically  used  in  compressible  fluid 
computations.  Essentially,  the  procedure  recognizes  that  for  a given 
polynomial  approximation  of  the  solution  of  degree  p over  an  element  of 
extent  h,  the  discretization  error  is  generally  0(hp+l).  At  certain  points 
within  the  element,  that  correspond  to  the  roots  of  the  difference  of  two 
Legendre  polynomials,  the  error  is  0(hp+2).  For  details  the  reader  is 
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directed  to  [1].  The  different  spatial  orders  of  accuracy  permit  a rapid 
a posteriori  estimate  of  the  accuracy  achieved  on  the  cells  at  any  one 
particular  grid  level.  Rather  than  obtaining  the  error  estimate  by  com- 
parison to  a test  integrating  in  time  on  a coarser  grid,  only  the  results 
from  a trial  time  step  on  the  current  grid  level  is  required. 

Changing  the  mathematical  error  estimation  method  is  not  enough. 
Tests  using  the  a posteriori  error  estimate  again  led  to  excessive  refine- 
ment and  consequent  loss  of  economy  of  computation.  An  examination 
of  the  flow  field  showed  that  this  occurred  even  in  regions  where  the 
fluid  turbulence  was  essentially  isotropic  and  suitable  for  modeling  by 
the  subgrid  scale  terms.  A physical  correction  of  the  error  estimator 
was  therefore  added.  In  the  regions  flagged  for  refinement  by  the  error 
estimator  two  physical  indicators  were  computed:  the  enstrophy  |V  x u\ 
and  the  scalar  product  of  the  local  velocity  and  vorticity  | ( V x u)  ■ u| . 
Only  when  one  of  the  physical  indicators  exceeded  a threshold  value 
was  the  region  subject  to  refinement.  In  effect  the  physical  indicators 
act  to  discriminate  interesting  dynamics  that  includes  prominent  vortex 
tubes  from  that  of  locally  near- isotropic  turbulence.  The  cutoff  value 
was  determined  by  numerical  experimentation.  This  is  unsatisfactory 
and  further  analysis  is  underway  to  establish  a more  rational  procedure 
of  establishing  a cutoff  value.  Other  indicators  suggested  by  coherent 
structure  eduction  procedures  are  also  undergoing  tests. 

1.5  Communication  between  grids 

An  important  aspect  in  a multilevel  algorithm  is  to  establish  proce- 
dures for  communication  of  data  between  grids  on  various  levels.  This 
involves  two  operations:  (1)  a prolongation  operation  P from  coarse  to 
fine  grids  and;  (2)  a restriction  operation  R from  fine  to  coarse  grids. 
Standard  AMR  typically  uses  a linear  interpolation  to  define  P and  a 
cell  averaging  procedure  to  define  R.  A modification  of  these  procedures 
was  found  to  be  necessary  in  the  context  of  applying  AMR  to  LES.  The 
prolongation  operator  is  applied  whenever  new  fine  grids  are  generated. 
Generally  there  is  significant  overlap  between  the  new  fine  grids  and 
previously  generated  grids  at  the  same  level  since  coherent  structures 
that  require  better  resolution  are  advected  with  the  mean  flow.  In  order 
not  to  lose  dynamic  content  (high  frequency  contributions),  fine  grids 
are  first  initialized  to  previously  computed  values  at  the  same  grid  level 
in  areas  of  grid  overlap.  The  prolongation  operator  is  only  applied  to 
newly  created  fine  grids,  at  level  l say,  where  no  grids  of  level  l existed 
previously.  Straightforward  linear  interpolation  was  found  to  induce 
excessive  subgrid  cell  excitation.  This  had  the  effect  of  increasing  the 
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number  iterations  in  the  pseudo-time  variable  (see  below)  required  to 
achieve  convergence  at  a given  physical  time  step.  Faster  convergence 
was  achieved  by  using  the  truncation  of  the  dispersion  relation  from  the 
Navier-Stokes  equations  to  the  wavenumbers  resolvable  on  the  newly  cre- 
ated grid.  This  relation  was  used  to  define  the  newly  resolvable  spectral 
content  by  assuming  that  the  energy  spectrum  follows  a power  law.  The 
power  law  coefficient  was  taken  to  depend  on  the  region  of  the  spectrum, 

1. e.  a quadratic  interpolation  of  the  power  law  coefficient  from  m — 4 in 
the  energy  containing  eddy  region  to  m — —5/3  in  the  inertial  range. 

A number  of  previous  multi-level  turbulence  simulations  (e.g.  [15]) 
mentioned  difficulties  at  the  boundaries  between  coarse  and  fine  grids 
due  to  the  effect  of  averaging.  Simple  averaging  was  found  to  have  the 
same  effect  in  the  computations  carried  out  here.  However,  conservative 
fix-ups  [5],  [4]  in  which  the  more  accurately  computed  fluxes  on  fine  grids 
are  used  to  update  adjoining  coarse  grid  values  was  found  to  work  well. 

2.  Numerical  Methods 

2.1  Unstructured  Grid  Algorithm 

The  unstructured  grid  algorithm  uses  an  embedded  tetrahedral  grid 
approach  implemented  as  an  octal  tree  structure  (OCTLES  - Octal  Tree 
Large  Eddy  Simulation  code)  [13].  There  are  several  options  for  treating 
the  convective  terms  from  the  Navier-Stokes  equations.  These  include 
multi-dimensional  fluctuation  splitting  [6],  second  order  reconstruction 
and  approximate  Riemann  solvers  along  the  cell  interface  normal  di- 
rection [17]  and  multi-dimensional  wave  propagation  techniques  [12]. 
The  convective  terms  are  advanced  in  time  explicitly.  Diffusive  terms 
are  treated  implicitly.  A pseudo-time  stepping  technique  [3]  is  used 
within  each  physical  time  step  to  solve  the  resulting  implicit  system. 
Subgrid-scale  turbulent  stresses  are  included  using  the  dynamic  model 
[10].  An  interesting  feature  of  this  application  of  the  dynamic  model  is 
that  subgrid-scale  stresses  are  computed  for  each  grid  pair  between  the 
coarsest,  initial  grid  and  the  finest  grid  present.  This  has  been  observed 
to  speed  up  the  convergence  of  the  pseudo-time  iterations  required  to 
determine  the  contribution  of  the  viscous  terms.  The  code  has  been  vali- 
dated [14]  by  comparison  to  a number  of  test  cases  proposed  in  [2],  Grid 
adaptation  is  carried  out  by  recursive  subdivision  of  an  initial  tetrahedral 
grid. 
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2.2  Structured  Cartesian  grid  algorithm 

An  alternative  to  the  complicated  program  structures  required  for 
unstructured  girds  is  to  use  embedded  Cartesian  grids  [5].  The  grid 
generation  procedures  are  much  simpler  in  this  case,  but  the  prob- 
lem of  boundary  representation  appears  for  general  geometries.  In  or- 
der to  test  the  AMR-LES  procedures,  sample  computations  have  been 
performed  only  in  a rectangular  geometry.  The  computer  code  for 
this  approach  (BEARCLAW  -Boundary  Embedded  Adaptive  Refine- 
ment, Conservation  Law)  allows  for  a number  of  treatments  of  the  con- 
vective and  source  terms.  These  parallel  those  presented  in  the  unstruc- 
tured grid  case.  The  BEARCLAW  code  may  be  freely  downloaded  from 
www.washington.edu/~claw. 

3.  Applications 

3.1  Stationary  and  rotational  channel  flow 

The  above  procedures  are  now  tested  on  channel  flow.  Both  a sta- 
tionary channel  and  one  rotating  along  an  axis  in  the  spanwise  mean 
flow  direction  are  considered.  The  AMR-LES  results  are  compared  to 
DNS  results  [11].  The  half  channel  width,  bulk  velocity  Reynolds  num- 
ber is  Re  = Ubh/u  — 2900.  The  rotating  channel  has  a Rossby  number 
Ro  = 2Clh/Ub  = 0.5.  The  second  order  wave  propagation  algorithms 
of  [12]  are  used  to  treat  the  convective  terms  in  both  the  structured 
and  unstructured  computations  presented  here.  A comparison  between 
the  DNS  results  and  those  obtained  by  the  two  AMR  methods  is  pre- 
sented in  fig.  1-2  for  the  turbulent  stress.  An  initial,  coarsest  level,  grid 
of  dimensions  32x32x32  was  used.  A tetrahedral  grid  was  obtained  by 
splitting  each  hexahedron  of  the  regular  Cartesian  grid  in  six  tetrahedra. 
The  AMR  results  are  within  1%  of  the  DNS  results  for  the  stationary 
channel  flow.  Both  the  structured  and  unstructured  methods  perform 
similarly  on  this  case.  This  is  thought  to  result  from  the  overall  sym- 
metry of  the  flow  and  the  fact  that  the  tetrahedral  grid  was  obtained 
by  destructuring  the  Cartesian  grid.  It  is  apparent  from  fig.  1 that  the 
AMR  algorithms  place  finer  resolution  in  the  high  shear  regions  closer 
to  the  walls.  When  rotation  is  applied  the  AMR-LES  methods  exhibit 
lower  accuracy.  The  maximum  error  observed  in  the  turbulent  stress 
was  2.4%.  Also,  the  loss  of  symmetry  imposed  by  the  rotation  led  to  dif- 
ferent grid  placement  in  the  unstructured  algorithm  with  respect  to  the 
Cartesian  algorithm.  Higher  error  levels  were  observed  for  the  tetrahe- 
dral grid  algorithm.  Instantaneous  vorticity  plots  of  the  flow  computed 
by  the  Cartesian  AMR-LES  method  are  presented  in  fig.  3-4.  It  may 
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be  observed  that  the  AMR  procedure  captures  the  rich  dynamics  of  the 
near-wall  region  and  the  relaminarization  which  occurs  on  the  rotating 
channel  suction  side. 

4.  Conclusions 

An  adaptive  mesh  refinement  methodology  for  large  eddy  simulation 
has  been  proposed.  LES  in  itself  reduces  the  dynamic  complexity  of 
simulating  fluid  turbulence.  It  is  asserted  that  adaptive  refinement  may 
offer  further  reductions  of  the  degrees  of  freedom  that  need  to  be  re- 
solved. The  AMR-LES  algorithm  proposed  here  differs  from  standard 
AMR  algorithms.  There  are  always  some  dynamics  taking  place  at  the 
unresolved  scale,  so  the  usual  error  estimation  procedures  used  to  add 
mesh  points  do  not  carry  over  directly.  Rather,  a combined  strategy 
incorporating  both  numerical  estimates  and  physical  flow  features  is  ap- 
plied. The  overall  effect  is  to  employ  greater  grid  resolution  in  regions 
in  which  the  subgrid  effects  relative  to  grid  level  l are  diagnosed  as  in- 
volving significant  deviations  from  isotropic  turbulence.  In  regions  in 
which  the  indicator  suggests  that  subgrid  fluctuations  are  isotropic  in 
nature,  a standard  dynamic,  eddy-viscosity  model  is  used  to  provide 
closure  terms.  Pull  testing  of  the  method  is  still  in  progress  and  shall  be 
reported  at  a later  date. 
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Figure  1.  Comparison  between  DNS  computation  (line)  [11],  Cartesian  grid  AMR 
(squares)  and  tetrahedral  grid  AMR  (diamonds).  The  turbulent  stress  along  the 
channel  span.  Notice  the  larger  grid  spacing  close  the  channel  centerline  and  the 
smaller  grid  spacing  near  the  walls. 
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Figure  2.  Similar  to  previous  figure  but  for  rotating  channel,  Ro  = 0.5. 


Figure  4-  Isovorticity  surfaces  computed  using  Cartesian  AMR,  rotating  channel 
Ro  = 0.5. 


